## selin code (original)
rm(list=ls())
library(MCMCpack)
library(rms)
library(plyr)
library(FactoMineR)
library(factoextra)
library(texreg)
library(ggmcmc)
library(systemfit)
library(AER)
library(stargazer)
library(lmtest)
library(foreign)

setwd("~/Dropbox/Projects/Patronage and Agency Independence")
selindata <- read.csv("Selin Agency Data.csv")

set.seed(1)
orig.estimates <- MCMCfactanal(~EOP+Cabinet+StatMandate+StatPermit+QuorumRules+Agency.specific.personnel+Term.Length+For.Cause+ServePresident+NumberMembers+Expertise+Conflict.of.Interest+Party.Balancing+Staggered.Terms+PAS.Head+President.Selects.Chair+Bureau+No.OMB.Budget.Review+No.OMB.Rule.Review+No.OMB.Communication.Review+Independent.Litigating+Independent.Funding+IG+Advisory.Commissions+Outside.Approval+Adjudication+ALJs,
                                data=selindata,
                                factors=2,
                                lambda.constraints=list(EOP=list(1,"-"),Cabinet=list(1,0),StatMandate=list(1,0),StatPermit=list(1,0),QuorumRules=list(1,0),Agency.specific.personnel=list(1,0),Term.Length=list(1,"+"),For.Cause=list(1,"+"),ServePresident=list(1,"-"),NumberMembers=list(1,0),Expertise=list(1,0),Conflict.of.Interest=list(1,0),Party.Balancing=list(1,0),Staggered.Terms=list(1,0),PAS.Head=list(1,0),President.Selects.Chair=list(1,0),Bureau=list(1,0),No.OMB.Budget.Review=list(2,0),No.OMB.Rule.Review=list(2,"+"),No.OMB.Communication.Review=list(2,"+"),Independent.Litigating=list(2,0),Independent.Funding=list(2,"+"),IG=list(2,0),Advisory.Commissions=list(2,0),Outside.Approval=list(2,0),Adjudication=list(2,0),ALJs=list(2,0)),
                                std.mean=TRUE, std.var=TRUE,verbose=100000,
                                mcmc = 1000000, burnin=100000,thin=1000,store.scores=TRUE)
saveRDS(orig.estimates, file = "orig-selin-estimates.rds")

# 
# 
# # next, estimate without coefficients from limitations on appointments (expertise/partisan balancing/conflict of interest)
new.estimates <- MCMCfactanal(~EOP+Cabinet+StatMandate+StatPermit+QuorumRules+Agency.specific.personnel+Term.Length+For.Cause+ServePresident+NumberMembers+Staggered.Terms+PAS.Head+President.Selects.Chair+Bureau+No.OMB.Budget.Review+No.OMB.Rule.Review+No.OMB.Communication.Review+Independent.Litigating+Independent.Funding+IG+Advisory.Commissions+Outside.Approval+Adjudication+ALJs,
                               data=selindata,
                               factors=2,
                               lambda.constraints=list(EOP=list(1,"-"),Cabinet=list(1,0),StatMandate=list(1,0),StatPermit=list(1,0),QuorumRules=list(1,0),Agency.specific.personnel=list(1,0),Term.Length=list(1,"+"),For.Cause=list(1,"+"),ServePresident=list(1,"-"),NumberMembers=list(1,0),Staggered.Terms=list(1,0),PAS.Head=list(1,0),President.Selects.Chair=list(1,0),Bureau=list(1,0),No.OMB.Budget.Review=list(2,0),No.OMB.Rule.Review=list(2,"+"),No.OMB.Communication.Review=list(2,"+"),Independent.Litigating=list(2,0),Independent.Funding=list(2,"+"),IG=list(2,0),Advisory.Commissions=list(2,0),Outside.Approval=list(2,0),Adjudication=list(2,0),ALJs=list(2,0)),
                               std.mean=TRUE, std.var=TRUE,verbose=100000,
                               mcmc = 1000000, burnin=100000,thin=1000,store.scores=TRUE)

saveRDS(new.estimates, file = "new-selin-estimates.rds")

exp.estimates <- MCMCfactanal(~EOP+Cabinet+StatMandate+StatPermit+QuorumRules+Agency.specific.personnel+Term.Length+For.Cause+ServePresident+NumberMembers+Conflict.of.Interest+Party.Balancing+Staggered.Terms+PAS.Head+President.Selects.Chair+Bureau+No.OMB.Budget.Review+No.OMB.Rule.Review+No.OMB.Communication.Review+Independent.Litigating+Independent.Funding+IG+Advisory.Commissions+Outside.Approval+Adjudication+ALJs,
                               data=selindata,
                               factors=2,
                               lambda.constraints=list(EOP=list(1,"-"),Cabinet=list(1,0),StatMandate=list(1,0),StatPermit=list(1,0),QuorumRules=list(1,0),Agency.specific.personnel=list(1,0),Term.Length=list(1,"+"),For.Cause=list(1,"+"),ServePresident=list(1,"-"),NumberMembers=list(1,0),Conflict.of.Interest=list(1,0),Party.Balancing=list(1,0),Staggered.Terms=list(1,0),PAS.Head=list(1,0),President.Selects.Chair=list(1,0),Bureau=list(1,0),No.OMB.Budget.Review=list(2,0),No.OMB.Rule.Review=list(2,"+"),No.OMB.Communication.Review=list(2,"+"),Independent.Litigating=list(2,0),Independent.Funding=list(2,"+"),IG=list(2,0),Advisory.Commissions=list(2,0),Outside.Approval=list(2,0),Adjudication=list(2,0),ALJs=list(2,0)),
                               std.mean=TRUE, std.var=TRUE,verbose=100000,
                               mcmc = 1000000, burnin=100000,thin=1000,store.scores=TRUE)

saveRDS(exp.estimates, file = "noexp-selin-estimates.rds")



orig.estimates <- readRDS(file = "orig-selin-estimates.rds")
new.estimates <- readRDS(file = "new-selin-estimates.rds")
exp.estimates <- readRDS(file = "noexp-selin-estimates.rds")


substrRight <- function(x, n){
  substr(x, nchar(x)-n+1, nchar(x))
}

phis.only.orig <- grep("phi", rownames(summary(orig.estimates)[[1]]))
phis.only <- grep("phi", rownames(summary(new.estimates)[[1]]))
phis.only.exp <- grep("phi", rownames(summary(exp.estimates)[[1]]))


agency.scores.orig <- summary(orig.estimates)[[1]][phis.only.orig,]
agency.scores <- summary(new.estimates)[[1]][phis.only,]
agency.scores.exp <- summary(exp.estimates)[[1]][phis.only.exp,]

factor.one <- agency.scores[which(substrRight(rownames(agency.scores),2) == "_1"),c(1,2)]
factor.two <- agency.scores[which(substrRight(rownames(agency.scores),2) == "_2"),c(1,2)]
factor.one.orig <- agency.scores.orig[which(substrRight(rownames(agency.scores.orig),2) == "_1"),c(1,2)]
factor.two.orig <- agency.scores.orig[which(substrRight(rownames(agency.scores.orig),2) == "_2"),c(1,2)]
factor.one.exp <- agency.scores.exp[which(substrRight(rownames(agency.scores.exp),2) == "_1"),c(1,2)]
factor.two.exp <- agency.scores.exp[which(substrRight(rownames(agency.scores.exp),2) == "_2"),c(1,2)]


new.selin <- data.frame(selindata, factor.one, factor.two, 
                        factor.one.orig, factor.two.orig, 
                        factor.one.exp, factor.two.exp, 
                        row.names = 1:345)
colnames(new.selin)[52:63] <- c("factor1_mean", "factor1_sd", "factor2_mean", "factor2_sd",
                                "factor1_mean_orig", "factor1_sd_orig", 
                                "factor2_mean_orig", "factor2_sd_orig",
                                "factor1_mean_exp", "factor1_sd_exp",
                                "factor2_mean_exp", "factor2_mean_exp") # factor 2 is independence
colnames(new.selin)[1] <- "FSAgency"


### scale the agency-level data from Hollibaugh-Horton-Lewis (AJPS)
agency.obs.temp <- read.csv("HHL Agency Data.csv")

agency.obs.temp <- subset(agency.obs.temp, select = c(fscode, fsagency, commission, 
                                                      numberofappointees, clintonlewisideo,
                                                      sotupriority, workforcesize, lnworkforcesize,
                                                      clerical, clericalpct, bluecollar,
                                                      bluecollarpct, clericalbluecollarpct,
                                                      professional, professionalpct))
colnames(agency.obs.temp)[4] <- "numappointees"
agency.obs.temp$lnprobias <- log(1 + agency.obs.temp$professionalpct) - 
                             log(1 + agency.obs.temp$clericalbluecollarpct)

appointee.obs <- read.csv("HHL Appointee Data.csv")
appointee.obs$PoliticsNew <- appointee.obs$PoliticsLastJob
appointee.obs$PoliticsNew[is.na(appointee.obs$PoliticsNew)] <- 0
appointee.obs$Postgrad <- appointee.obs$MD.Mphil + appointee.obs$Masters + appointee.obs$PhD

appointee.obs$Department[grep("Export-Import Bank", appointee.obs$Department)] <- "Export-Import Bank"

agency.obs2 <- ddply(appointee.obs, .variables = .(fscode),
                        summarise,
                        numappointees = mean(NumberofAppointeesinAgency, na.rm = TRUE),
                        avgcampaign = mean(Campaign, na.rm = TRUE),
                        avgtransition = mean(Transition, na.rm = TRUE),
                        avgcamptrans = mean(I(Campaign == 1|Transition == 1), 
                                                         na.rm = TRUE),
                        avgbush = mean(Bush, na.rm = TRUE),
                        avgclinton = mean(Clinton, na.rm = TRUE),
                        avgbushorclinton = mean(I(Bush == 1|Clinton == 1), 
                                                na.rm = TRUE),
                        avgagency = mean(Agency, na.rm = TRUE),
                        avgpostgrad = mean(Postgrad, na.rm = TRUE),
                        avgbachelors = mean(Bachelors, na.rm = TRUE),
                        avgphd = mean(PhD, na.rm = TRUE),
                        avgsubjectall = mean(SubjectAll, na.rm = TRUE),
                        avgpolitics = mean(PoliticsNew, na.rm = TRUE),
                        avgconnection = mean(Connection, na.rm = TRUE),
                        avgharvard = mean(Harvard, na.rm = TRUE),
                        avggovt = mean(Gov.t, na.rm = TRUE),
                        avgschedc = mean(SchedC, na.rm = TRUE))

agency.obs2 <- data.frame(subset(agency.obs2, numappointees != 0), row.names = NULL) 
avgschedc_raw <- agency.obs2[,dim(agency.obs2)[2]]
agency.obs2 <- cbind(agency.obs2[,1:2], apply(agency.obs2[,-c(1:2)],2,scale), avgschedc_raw)
patron.dim <- with(agency.obs2,
                   PCA(cbind(avgpolitics, avgcamptrans, avgconnection),row.w= numappointees,
                       graph = FALSE))   
expert.dim <- with(agency.obs2,
                   PCA(cbind(avgpostgrad, avgbachelors, avgbushorclinton, 
                             avgagency, avgsubjectall, avggovt), row.w = numappointees,
                       graph = FALSE))
both.dim <- with(agency.obs2,
                 PCA(cbind(avgpostgrad, avgbachelors, avgbushorclinton, 
                           avgagency, avgsubjectall, avggovt,
                           avgpolitics, avgcamptrans, avgconnection), 
                     row.w = numappointees,
                     graph = FALSE))

agency.obs <- data.frame(agency.obs2, 
                         latentpatron = scale(patron.dim$ind$coord[,1]),
                         latentexpert = scale(expert.dim$ind$coord[,1]),
                         patronexpert = -scale(both.dim$ind$coord[,1]))                        
colnames(agency.obs.temp)[2] <- "FSAgency"

agency.obs <- merge(agency.obs.temp, agency.obs, by = "fscode")                        

agency.walkthrough <- subset(agency.obs, select = c(fscode, FSAgency))
colnames(agency.walkthrough)[1] <- "fscode"
appointee.obs <- merge(agency.walkthrough, appointee.obs)

colnames(appointee.obs)[2] <- "FSAgency"
colnames(agency.obs)[2] <- "FSAgency"
colnames(agency.obs)[1] <- "fscode"

new.selin$FSAgency <- as.character(new.selin$FSAgency)
new.selin$FSAgency[which(new.selin$FSAgency == "Export-Import Bank of the United States")] <- "Export-Import Bank"
new.selin$FSAgency[which(new.selin$FSAgency == "Commodities Futures Trading Commission")] <- "Commodity Futures Trading Commission"
new.selin$FSAgency[which(new.selin$FSAgency == "United States Trade and Development Agency")] <- "Trade and Development Agency"
new.selin$FSAgency[which(new.selin$FSAgency == "Office of Science and Technology")] <- "Office of Science and Technology Policy"


selin.agency.merged <- merge(new.selin, agency.obs)
selin.appointee.merged <- merge(new.selin, appointee.obs, by = "FSAgency")
selin.agency.merged$factor1_mean <- scale(selin.agency.merged$factor1_mean)
selin.agency.merged$factor2_mean <- scale(selin.agency.merged$factor2_mean)
selin.agency.merged$factor1_mean_orig <- scale(selin.agency.merged$factor1_mean_orig)
selin.agency.merged$factor2_mean_orig <- scale(selin.agency.merged$factor2_mean_orig)
selin.agency.merged$factor1_mean_exp <- scale(selin.agency.merged$factor1_mean_exp)
selin.agency.merged$factor2_mean_exp <- scale(selin.agency.merged$factor2_mean_exp)
selin.agency.merged$latentpatron <- scale(selin.agency.merged$latentpatron)
selin.agency.merged$latentexpert <- scale(selin.agency.merged$latentexpert)
selin.agency.merged$patronexpert <- scale(selin.agency.merged$patronexpert)

### factor one is policy independence
### factor two is decisionmaker independence



### Revised Analysis
selin.agency.merged$bluecollarpct_scale <- scale(selin.agency.merged$bluecollarpct)
selin.agency.merged$clericalpct_scale <- scale(selin.agency.merged$clericalpct)
selin.agency.merged$professional_scale <- scale(selin.agency.merged$professionalpct)
selin.agency.merged$numappointees <- selin.agency.merged$numappointees.y

set.seed(1)
mod.both.1 <- ols(latentpatron ~ factor1_mean + factor2_mean, 
                  data = selin.agency.merged, 
                  x = TRUE, y = TRUE)
mod.both.1.robust <- robcov(mod.both.1)
mod.both.1.boot <- bootcov(mod.both.1, B = 10000)


mod.both.1.orig <- ols(latentpatron ~ factor1_mean_orig + factor2_mean_orig, 
                       data = selin.agency.merged, 
                       x = TRUE, y = TRUE)
mod.both.1.orig.robust <- robcov(mod.both.1.orig)
mod.both.1.orig.boot <- bootcov(mod.both.1.orig, B = 10000)

mod.both.1.exp <- ols(latentpatron ~ factor1_mean_exp + factor2_mean_exp, 
                       data = selin.agency.merged, 
                       x = TRUE, y = TRUE)
mod.both.1.exp.robust <- robcov(mod.both.1.exp)
mod.both.1.exp.boot <- bootcov(mod.both.1.exp, B = 10000)


mod.both.1.weight <- ols(latentpatron ~ factor1_mean + factor2_mean, 
                         data = selin.agency.merged, x = TRUE, y = TRUE,
                         weights = numappointees)
mod.both.1.weight.robust <- robcov(mod.both.1.weight)
mod.both.1.weight.boot <- bootcov(mod.both.1.weight, B = 10000)

mod.both.1.orig.weight <- ols(latentpatron ~ factor1_mean_orig + factor2_mean_orig, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.both.1.orig.weight.robust <- robcov(mod.both.1.orig.weight)
mod.both.1.orig.weight.boot <- bootcov(mod.both.1.orig.weight, B = 10000)

mod.both.1.exp.weight <- ols(latentpatron ~ factor1_mean_exp + factor2_mean_exp, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.both.1.exp.weight.robust <- robcov(mod.both.1.exp.weight)
mod.both.1.exp.weight.boot <- bootcov(mod.both.1.exp.weight, B = 10000)


mod.all.1 <- ols(latentpatron ~ factor1_mean + factor2_mean + 
                   clintonlewisideo + lnworkforcesize + lnprobias + 
                   sotupriority, 
                 data = selin.agency.merged, x = TRUE, y = TRUE)
mod.all.1.robust <- robcov(mod.all.1)
mod.all.1.boot <- bootcov(mod.all.1, B = 10000)


mod.all.1.orig <- ols(latentpatron ~ factor1_mean_orig + factor2_mean_orig + 
                        clintonlewisideo + lnworkforcesize + lnprobias + 
                        sotupriority, 
                      data = selin.agency.merged, x = TRUE, y = TRUE)
mod.all.1.orig.robust <- robcov(mod.all.1.orig)
mod.all.1.orig.boot <- bootcov(mod.all.1.orig, B = 10000)

mod.all.1.exp <- ols(latentpatron ~ factor1_mean_exp + factor2_mean_exp + 
                        clintonlewisideo + lnworkforcesize + lnprobias + 
                        sotupriority, 
                      data = selin.agency.merged, x = TRUE, y = TRUE)
mod.all.1.exp.robust <- robcov(mod.all.1.exp)
mod.all.1.exp.boot <- bootcov(mod.all.1.exp, B = 10000)


mod.all.1.weight <- ols(latentpatron ~ factor1_mean + factor2_mean + 
                          clintonlewisideo + lnworkforcesize + lnprobias + 
                          sotupriority, 
                        data = selin.agency.merged, x = TRUE, y = TRUE,
                        weights = numappointees)
mod.all.1.weight.robust <- robcov(mod.all.1.weight)
mod.all.1.weight.boot <- bootcov(mod.all.1.weight, B = 10000)

mod.all.1.orig.weight <- ols(latentpatron ~ factor1_mean_orig + factor2_mean_orig + 
                               clintonlewisideo + lnworkforcesize + lnprobias + 
                               sotupriority, 
                             data = selin.agency.merged, x = TRUE, y = TRUE,
                             weights = numappointees)
mod.all.1.orig.weight.robust <- robcov(mod.all.1.orig.weight)
mod.all.1.orig.weight.boot <- bootcov(mod.all.1.orig.weight, B = 10000)

mod.all.1.exp.weight <- ols(latentpatron ~ factor1_mean_exp + factor2_mean_exp + 
                               clintonlewisideo + lnworkforcesize + lnprobias + 
                               sotupriority, 
                             data = selin.agency.merged, x = TRUE, y = TRUE,
                             weights = numappointees)
mod.all.1.exp.weight.robust <- robcov(mod.all.1.exp.weight)
mod.all.1.exp.weight.boot <- bootcov(mod.all.1.exp.weight, B = 10000)


mod.uber.1 <- ols(latentpatron ~ factor1_mean + factor2_mean + 
                    clintonlewisideo + lnworkforcesize + 
                    sotupriority + bluecollarpct_scale + 
                    clericalpct_scale + professional_scale, 
                  data = selin.agency.merged, x = TRUE, y = TRUE)
mod.uber.1.robust <- robcov(mod.uber.1)
mod.uber.1.boot <- bootcov(mod.uber.1, B = 10000)

mod.uber.1.orig <- ols(latentpatron ~ factor1_mean_orig + factor2_mean_orig + 
                         clintonlewisideo + lnworkforcesize + 
                         sotupriority + bluecollarpct_scale + 
                         clericalpct_scale + professional_scale, 
                       data = selin.agency.merged, x = TRUE, y = TRUE)
mod.uber.1.orig.robust <- robcov(mod.uber.1.orig)
mod.uber.1.orig.boot <- bootcov(mod.uber.1.orig, B = 10000)

mod.uber.1.exp <- ols(latentpatron ~ factor1_mean_exp + factor2_mean_exp + 
                         clintonlewisideo + lnworkforcesize + 
                         sotupriority + bluecollarpct_scale + 
                         clericalpct_scale + professional_scale, 
                       data = selin.agency.merged, x = TRUE, y = TRUE)
mod.uber.1.exp.robust <- robcov(mod.uber.1.exp)
mod.uber.1.exp.boot <- bootcov(mod.uber.1.exp, B = 10000)


mod.uber.1.weight <- ols(latentpatron ~ factor1_mean + factor2_mean + 
                           clintonlewisideo + lnworkforcesize + 
                           sotupriority + bluecollarpct_scale + 
                           clericalpct_scale + professional_scale, 
                         data = selin.agency.merged, x = TRUE, y = TRUE,
                         weights = numappointees)
mod.uber.1.weight.robust <- robcov(mod.uber.1.weight)
mod.uber.1.weight.boot <- bootcov(mod.uber.1.weight, B = 10000)

mod.uber.1.orig.weight <- ols(latentpatron ~ factor1_mean_orig + factor2_mean_orig + 
                                clintonlewisideo + lnworkforcesize + 
                                sotupriority + bluecollarpct_scale + 
                                clericalpct_scale + professional_scale, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.uber.1.orig.weight.robust <- robcov(mod.uber.1.orig.weight)
mod.uber.1.orig.weight.boot <- bootcov(mod.uber.1.orig.weight, B = 10000)

mod.uber.1.exp.weight <- ols(latentpatron ~ factor1_mean_exp + factor2_mean_exp + 
                                clintonlewisideo + lnworkforcesize + 
                                sotupriority + bluecollarpct_scale + 
                                clericalpct_scale + professional_scale, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.uber.1.exp.weight.robust <- robcov(mod.uber.1.exp.weight)
mod.uber.1.exp.weight.boot <- bootcov(mod.uber.1.exp.weight, B = 10000)


mod.both.2 <- ols(patronexpert ~ factor1_mean + factor2_mean, 
                  data = selin.agency.merged, 
                  x = TRUE, y = TRUE)
mod.both.2.robust <- robcov(mod.both.2)
mod.both.2.boot <- bootcov(mod.both.2, B = 10000)


mod.both.2.orig <- ols(patronexpert ~ factor1_mean_orig + factor2_mean_orig, 
                       data = selin.agency.merged, 
                       x = TRUE, y = TRUE)
mod.both.2.orig.robust <- robcov(mod.both.2.orig)
mod.both.2.orig.boot <- bootcov(mod.both.2.orig, B = 10000)

mod.both.2.exp <- ols(patronexpert ~ factor1_mean_exp + factor2_mean_exp, 
                       data = selin.agency.merged, 
                       x = TRUE, y = TRUE)
mod.both.2.exp.robust <- robcov(mod.both.2.exp)
mod.both.2.exp.boot <- bootcov(mod.both.2.exp, B = 10000)


mod.both.2.weight <- ols(patronexpert ~ factor1_mean + factor2_mean, 
                         data = selin.agency.merged, x = TRUE, y = TRUE,
                         weights = numappointees)
mod.both.2.weight.robust <- robcov(mod.both.2.weight)
mod.both.2.weight.boot <- bootcov(mod.both.2.weight, B = 10000)

mod.both.2.orig.weight <- ols(patronexpert ~ factor1_mean_orig + factor2_mean_orig, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.both.2.orig.weight.robust <- robcov(mod.both.2.orig.weight)
mod.both.2.orig.weight.boot <- bootcov(mod.both.2.orig.weight, B = 10000)

mod.both.2.exp.weight <- ols(patronexpert ~ factor1_mean_exp + factor2_mean_exp, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.both.2.exp.weight.robust <- robcov(mod.both.2.exp.weight)
mod.both.2.exp.weight.boot <- bootcov(mod.both.2.exp.weight, B = 10000)


mod.all.2 <- ols(patronexpert ~ factor1_mean + factor2_mean + 
                   clintonlewisideo + lnworkforcesize + lnprobias + 
                   sotupriority, 
                 data = selin.agency.merged, x = TRUE, y = TRUE)
mod.all.2.robust <- robcov(mod.all.2)
mod.all.2.boot <- bootcov(mod.all.2, B = 10000)


mod.all.2.orig <- ols(patronexpert ~ factor1_mean_orig + factor2_mean_orig + 
                        clintonlewisideo + lnworkforcesize + lnprobias + 
                        sotupriority, 
                      data = selin.agency.merged, x = TRUE, y = TRUE)
mod.all.2.orig.robust <- robcov(mod.all.2.orig)
mod.all.2.orig.boot <- bootcov(mod.all.2.orig, B = 10000)

mod.all.2.exp <- ols(patronexpert ~ factor1_mean_exp + factor2_mean_exp + 
                        clintonlewisideo + lnworkforcesize + lnprobias + 
                        sotupriority, 
                      data = selin.agency.merged, x = TRUE, y = TRUE)
mod.all.2.exp.robust <- robcov(mod.all.2.exp)
mod.all.2.exp.boot <- bootcov(mod.all.2.exp, B = 10000)


mod.all.2.weight <- ols(patronexpert ~ factor1_mean + factor2_mean + 
                          clintonlewisideo + lnworkforcesize + lnprobias + 
                          sotupriority, 
                        data = selin.agency.merged, x = TRUE, y = TRUE,
                        weights = numappointees)
mod.all.2.weight.robust <- robcov(mod.all.2.weight)
mod.all.2.weight.boot <- bootcov(mod.all.2.weight, B = 10000)

mod.all.2.orig.weight <- ols(patronexpert ~ factor1_mean_orig + factor2_mean_orig + 
                               clintonlewisideo + lnworkforcesize + lnprobias + 
                               sotupriority, 
                             data = selin.agency.merged, x = TRUE, y = TRUE,
                             weights = numappointees)
mod.all.2.orig.weight.robust <- robcov(mod.all.2.orig.weight)
mod.all.2.orig.weight.boot <- bootcov(mod.all.2.orig.weight, B = 10000)

mod.all.2.exp.weight <- ols(patronexpert ~ factor1_mean_exp + factor2_mean_exp + 
                               clintonlewisideo + lnworkforcesize + lnprobias + 
                               sotupriority, 
                             data = selin.agency.merged, x = TRUE, y = TRUE,
                             weights = numappointees)
mod.all.2.exp.weight.robust <- robcov(mod.all.2.exp.weight)
mod.all.2.exp.weight.boot <- bootcov(mod.all.2.exp.weight, B = 10000)


mod.uber.2 <- ols(patronexpert ~ factor1_mean + factor2_mean + 
                    clintonlewisideo + lnworkforcesize + 
                    sotupriority + bluecollarpct_scale + 
                    clericalpct_scale + professional_scale, 
                  data = selin.agency.merged, x = TRUE, y = TRUE)
mod.uber.2.robust <- robcov(mod.uber.2)
mod.uber.2.boot <- bootcov(mod.uber.2, B = 10000)

mod.uber.2.orig <- ols(patronexpert ~ factor1_mean_orig + factor2_mean_orig + 
                         clintonlewisideo + lnworkforcesize + 
                         sotupriority + bluecollarpct_scale + 
                         clericalpct_scale + professional_scale, 
                       data = selin.agency.merged, x = TRUE, y = TRUE)
mod.uber.2.orig.robust <- robcov(mod.uber.2.orig)
mod.uber.2.orig.boot <- bootcov(mod.uber.2.orig, B = 10000)

mod.uber.2.exp <- ols(patronexpert ~ factor1_mean_exp + factor2_mean_exp + 
                         clintonlewisideo + lnworkforcesize + 
                         sotupriority + bluecollarpct_scale + 
                         clericalpct_scale + professional_scale, 
                       data = selin.agency.merged, x = TRUE, y = TRUE)
mod.uber.2.exp.robust <- robcov(mod.uber.2.exp)
mod.uber.2.exp.boot <- bootcov(mod.uber.2.exp, B = 10000)



mod.uber.2.weight <- ols(patronexpert ~ factor1_mean + factor2_mean + 
                           clintonlewisideo + lnworkforcesize + 
                           sotupriority + bluecollarpct_scale + 
                           clericalpct_scale + professional_scale, 
                         data = selin.agency.merged, x = TRUE, y = TRUE,
                         weights = numappointees)
mod.uber.2.weight.robust <- robcov(mod.uber.2.weight)
mod.uber.2.weight.boot <- bootcov(mod.uber.2.weight, B = 10000)

mod.uber.2.orig.weight <- ols(patronexpert ~ factor1_mean_orig + factor2_mean_orig + 
                                clintonlewisideo + lnworkforcesize + 
                                sotupriority + bluecollarpct_scale + 
                                clericalpct_scale + professional_scale, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.uber.2.orig.weight.robust <- robcov(mod.uber.2.orig.weight)
mod.uber.2.orig.weight.boot <- bootcov(mod.uber.2.orig.weight, B = 10000)


mod.uber.2.exp.weight <- ols(patronexpert ~ factor1_mean_exp + factor2_mean_exp + 
                                clintonlewisideo + lnworkforcesize + 
                                sotupriority + bluecollarpct_scale + 
                                clericalpct_scale + professional_scale, 
                              data = selin.agency.merged, x = TRUE, y = TRUE,
                              weights = numappointees)
mod.uber.2.exp.weight.robust <- robcov(mod.uber.2.exp.weight)
mod.uber.2.exp.weight.boot <- bootcov(mod.uber.2.exp.weight, B = 10000)




# OLS output function, for the nominee-president ideological divergence models.
extract.ols <- function(model, include.nobs = TRUE,
                        include.rsquared = TRUE, include.adjrs = TRUE,
                        include.lr = TRUE, ...){
  names <- dimnames(model$var)[[1]]
  co <- c(model$coefficients)
  se <- sqrt(diag(model$var))
  pval <- 2*(1-pnorm(abs(co/se)))
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.nobs == TRUE) {
    n <- nobs(model)
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  if (include.rsquared == TRUE) {
    rs <- model$stats["R2"]
    gof <- c(gof, rs)
    gof.names <- c(gof.names, "R$^2$")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.adjrs == TRUE) {
    adj <- 1 - (1 - model$stats["R2"]) * (n - 1) / (n - model$stats["d.f."] - 1)
    gof <- c(gof, adj)
    gof.names <- c(gof.names, "Adjusted R$^2$")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.lr == TRUE) {
    LR <- model$stats["Model L.R."]
    LRT <- as.numeric(pchisq(LR, model$stats['d.f.'], lower.tail = FALSE))
    gof <- c(gof, LR, LRT)
    gof.names <- c(gof.names, "Likelihood Ratio Test", "Likelihood Ratio Test p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  tr <- createTexreg(coef.names = names, coef = co, se = se, 
                     pvalues = pval, gof.names = gof.names, gof = gof, gof.decimal = gof.decimal)
  return(tr)
}


### using largest standard errors across standard/bootstrapped/sandwich robust
extract.big <- function(model){
  foo <- extract(eval(parse(text = paste(model))))
  foo@se <- apply(rbind(extract(eval(parse(text = paste(model, sep = ""))))@se, 
                        extract(eval(parse(text = paste(model,".robust", sep = ""))))@se,
                        extract(eval(parse(text = paste(model,".boot", sep = ""))))@se),
                  2,max)
  foo@pvalues <- apply(rbind(extract(eval(parse(text = paste(model, sep = ""))))@pvalues, 
                             extract(eval(parse(text = paste(model,".robust", sep = ""))))@pvalues,
                             extract(eval(parse(text = paste(model,".boot", sep = ""))))@pvalues),
                       2,max)
  return(foo)
}

mod.both.1.extract <- extract.big("mod.both.1")
mod.both.1.weight.extract <- extract.big("mod.both.1.weight")
mod.all.1.extract <- extract.big("mod.all.1")
mod.all.1.weight.extract <- extract.big("mod.all.1.weight")
mod.uber.1.extract <- extract.big("mod.uber.1")
mod.uber.1.weight.extract <- extract.big("mod.uber.1.weight")

mod.both.1.orig.extract <- extract.big("mod.both.1.orig")
mod.both.1.orig.weight.extract <- extract.big("mod.both.1.orig.weight")
mod.all.1.orig.extract <- extract.big("mod.all.1.orig")
mod.all.1.orig.weight.extract <- extract.big("mod.all.1.orig.weight")
mod.uber.1.orig.extract <- extract.big("mod.uber.1.orig")
mod.uber.1.orig.weight.extract <- extract.big("mod.uber.1.orig.weight")

mod.both.1.exp.extract <- extract.big("mod.both.1.exp")
mod.both.1.exp.weight.extract <- extract.big("mod.both.1.exp.weight")
mod.all.1.exp.extract <- extract.big("mod.all.1.exp")
mod.all.1.exp.weight.extract <- extract.big("mod.all.1.exp.weight")
mod.uber.1.exp.extract <- extract.big("mod.uber.1.exp")
mod.uber.1.exp.weight.extract <- extract.big("mod.uber.1.exp.weight")


mod.both.2.extract <- extract.big("mod.both.2")
mod.both.2.weight.extract <- extract.big("mod.both.2.weight")
mod.all.2.extract <- extract.big("mod.all.2")
mod.all.2.weight.extract <- extract.big("mod.all.2.weight")
mod.uber.2.extract <- extract.big("mod.uber.2")
mod.uber.2.weight.extract <- extract.big("mod.uber.2.weight")

mod.both.2.orig.extract <- extract.big("mod.both.2.orig")
mod.both.2.orig.weight.extract <- extract.big("mod.both.2.orig.weight")
mod.all.2.orig.extract <- extract.big("mod.all.2.orig")
mod.all.2.orig.weight.extract <- extract.big("mod.all.2.orig.weight")
mod.uber.2.orig.extract <- extract.big("mod.uber.2.orig")
mod.uber.2.orig.weight.extract <- extract.big("mod.uber.2.orig.weight")

mod.both.2.exp.extract <- extract.big("mod.both.2.exp")
mod.both.2.exp.weight.extract <- extract.big("mod.both.2.exp.weight")
mod.all.2.exp.extract <- extract.big("mod.all.2.exp")
mod.all.2.exp.weight.extract <- extract.big("mod.all.2.exp.weight")
mod.uber.2.exp.extract <- extract.big("mod.uber.2.exp")
mod.uber.2.exp.weight.extract <- extract.big("mod.uber.2.exp.weight")


texreg(list(mod.both.1.extract, mod.all.1.extract, mod.uber.1.extract,
            mod.both.1.weight.extract, mod.all.1.weight.extract, mod.uber.1.weight.extract),
       digits = 3,
       custom.coef.names = c("Constant",
                             "Policy Independence",
                             "Decisionmaker Independence",
                             "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency",
                             "Blue-Collar Employee Proportion",
                             "Clerical Employee Proportion",
                             "Professional Employee Proportion"),
       reorder.coef = c(3,2,4:10,1),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       custom.note = "",
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Model B-1",
                              "Model B-2",
                              "Model B-3",
                              "Model B-4",
                              "Model B-5",
                              "Model B-6"),
       file = "Table-B1.tex")

texreg(list(mod.both.1.orig.extract, 
            mod.all.1.orig.extract, 
            mod.uber.1.orig.extract,
            mod.both.1.orig.weight.extract, 
            mod.all.1.orig.weight.extract, 
            mod.uber.1.orig.weight.extract),
       digits = 3,
       custom.coef.names = c("Constant",
                             "Policy Independence",
                             "Decisionmaker Independence",
                             "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency",
                             "Blue-Collar Employee Proportion",
                             "Clerical Employee Proportion",
                             "Professional Employee Proportion"),
       reorder.coef = c(3,2,4:10,1),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       custom.note = "",
       stars = c(0.01, 0.05, 0.1), 
       file = "Table-1.tex")


texreg(list(mod.both.1.exp.extract, 
            mod.all.1.exp.extract, 
            mod.uber.1.exp.extract,
            mod.both.1.exp.weight.extract, 
            mod.all.1.exp.weight.extract, 
            mod.uber.1.exp.weight.extract),
       digits = 3,
       custom.coef.names = c("Constant",
                             "Policy Independence",
                             "Decisionmaker Independence",
                             "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency",
                             "Blue-Collar Employee Proportion",
                             "Clerical Employee Proportion",
                             "Professional Employee Proportion"),
       reorder.coef = c(3,2,4:10,1),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       custom.note = "",
       stars = c(0.01, 0.05, 0.1), 
       custom.model.names = c("Model C-1",
                              "Model C-2",
                              "Model C-3",
                              "Model C-4",
                              "Model C-5",
                              "Model C-6"),
       file = "Table-C1.tex")




texreg(list(mod.both.2.extract, mod.all.2.extract, mod.uber.2.extract,
            mod.both.2.weight.extract, mod.all.2.weight.extract, mod.uber.2.weight.extract),
       digits = 3,
       custom.coef.names = c("Constant",
                             "Policy Independence",
                             "Decisionmaker Independence",
                             "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency",
                             "Blue-Collar Employee Proportion",
                             "Clerical Employee Proportion",
                             "Professional Employee Proportion"),
       reorder.coef = c(3,2,4:10,1),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       custom.note = "",
       stars = c(0.01, 0.05, 0.1), 
       custom.model.names = c("Model E-7",
                              "Model E-8",
                              "Model E-9",
                              "Model E-10",
                              "Model E-11",
                              "Model E-12"),
       file = "Table-E2.tex")


texreg(list(mod.both.2.exp.extract, 
            mod.all.2.exp.extract, 
            mod.uber.2.exp.extract,
            mod.both.2.exp.weight.extract, 
            mod.all.2.exp.weight.extract, 
            mod.uber.2.exp.weight.extract),
       digits = 3,
       custom.coef.names = c("Constant",
                             "Policy Independence",
                             "Decisionmaker Independence",
                             "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency",
                             "Blue-Collar Employee Proportion",
                             "Clerical Employee Proportion",
                             "Professional Employee Proportion"),
       reorder.coef = c(3,2,4:10,1),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       custom.note = "",
       stars = c(0.01, 0.05, 0.1), 
       custom.model.names = c("Model E-13",
                              "Model E-14",
                              "Model E-15",
                              "Model E-16",
                              "Model E-17",
                              "Model E-18"),
       file = "Table-E3.tex")


texreg(list(mod.both.2.orig.extract, 
            mod.all.2.orig.extract, 
            mod.uber.2.orig.extract,
            mod.both.2.orig.weight.extract, 
            mod.all.2.orig.weight.extract, 
            mod.uber.2.orig.weight.extract),
       digits = 3,
       custom.coef.names = c("Constant",
                             "Policy Independence",
                             "Decisionmaker Independence",
                             "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency",
                             "Blue-Collar Employee Proportion",
                             "Clerical Employee Proportion",
                             "Professional Employee Proportion"),
       reorder.coef = c(3,2,4:10,1),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       custom.note = "",
       stars = c(0.01, 0.05, 0.1), 
       custom.model.names = c("Model E-1",
                              "Model E-2",
                              "Model E-3",
                              "Model E-4",
                              "Model E-5",
                              "Model E-6"),
       file = "Table-E1.tex")



### make 2x2 facet plot by rbinding everything together
selin.plot.frame <- rbind(data.frame(selin.agency.merged, 
                                     Type = "Decisionmaker Independence",
                                     Dimension =selin.agency.merged$factor2_mean,
                                     Weighted = "Unweighted",
                                     Weight = 1),
                          data.frame(selin.agency.merged, 
                                     Type = "Decisionmaker Independence",
                                     Dimension = selin.agency.merged$factor2_mean,
                                     Weighted = "Weighted",
                                     Weight = selin.agency.merged$numappointees.x),
                          data.frame(selin.agency.merged, 
                                     Type = "Policy Independence",
                                     Dimension = selin.agency.merged$factor1_mean,
                                     Weighted = "Unweighted",
                                     Weight = 1),
                          data.frame(selin.agency.merged, 
                                     Type = "Policy Independence",
                                     Dimension = selin.agency.merged$factor1_mean,
                                     Weighted = "Weighted",
                                     Weight = selin.agency.merged$numappointees.x))

selin.plot.frame.orig <- rbind(data.frame(selin.agency.merged, 
                                     Type = "Decisionmaker Independence",
                                     Dimension =selin.agency.merged$factor2_mean_orig,
                                     Weighted = "Unweighted",
                                     Weight = 1),
                          data.frame(selin.agency.merged, 
                                     Type = "Decisionmaker Independence",
                                     Dimension = selin.agency.merged$factor2_mean_orig,
                                     Weighted = "Weighted",
                                     Weight = selin.agency.merged$numappointees.x),
                          data.frame(selin.agency.merged, 
                                     Type = "Policy Independence",
                                     Dimension = selin.agency.merged$factor1_mean_orig,
                                     Weighted = "Unweighted",
                                     Weight = 1),
                          data.frame(selin.agency.merged, 
                                     Type = "Policy Independence",
                                     Dimension = selin.agency.merged$factor1_mean_orig,
                                     Weighted = "Weighted",
                                     Weight = selin.agency.merged$numappointees.x))


selin.plot.frame.exp <- rbind(data.frame(selin.agency.merged, 
                                         Type = "Decisionmaker Independence",
                                         Dimension =selin.agency.merged$factor2_mean_exp,
                                         Weighted = "Unweighted",
                                         Weight = 1),
                              data.frame(selin.agency.merged, 
                                         Type = "Decisionmaker Independence",
                                         Dimension = selin.agency.merged$factor2_mean_exp,
                                         Weighted = "Weighted",
                                         Weight = selin.agency.merged$numappointees.x),
                              data.frame(selin.agency.merged, 
                                         Type = "Policy Independence",
                                         Dimension = selin.agency.merged$factor1_mean_exp,
                                         Weighted = "Unweighted",
                                         Weight = 1),
                              data.frame(selin.agency.merged, 
                                         Type = "Policy Independence",
                                         Dimension = selin.agency.merged$factor1_mean_exp,
                                         Weighted = "Weighted",
                                         Weight = selin.agency.merged$numappointees.x))



g1.big <- ggplot(selin.plot.frame,  aes(x = Dimension, y = latentpatron, 
                                        weight = Weight)) + 
  geom_point(alpha = I(0.3)) + 
  stat_smooth(color = "black", level = 0.9, span = 1.2) + 
  theme_bw() + 
  xlab('\nIndependence\n(Number of Standard Deviations from Mean)\n') + 
  ylab('Patronage Dimension\n(Number of Standard Deviations from Mean)\n') + 
  facet_grid(Weighted~Type, scales = "free_x")

g1.big.both <- ggplot(selin.plot.frame,  aes(x = Dimension, y = patronexpert, 
                                             weight = Weight)) + 
  geom_point(alpha = I(0.3)) + 
  stat_smooth(color = "black", level = 0.9, span = 1.2) + 
  theme_bw() + 
  xlab('\nIndependence\n(Number of Standard Deviations from Mean)\n') + 
  ylab('Patronage-Expertise Dimension\n(Number of Standard Deviations from Mean)\n') + 
  facet_grid(Weighted~Type, scales = "free")



g1.big.orig <- ggplot(selin.plot.frame.orig,  aes(x = Dimension, y = latentpatron, 
                                        weight = Weight)) + 
  geom_point(alpha = I(0.3)) + 
  stat_smooth(color = "black", level = 0.9, span = 1.2) + 
  theme_bw() + 
  xlab('\nIndependence\n(Number of Standard Deviations from Mean)\n') + 
  ylab('Patronage Dimension\n(Number of Standard Deviations from Mean)\n') + 
  facet_grid(Weighted~Type, scales = "free_x")

g1.big.both.orig <- ggplot(selin.plot.frame.orig,  aes(x = Dimension, y = patronexpert, 
                                             weight = Weight)) + 
  geom_point(alpha = I(0.3)) + 
  stat_smooth(color = "black", level = 0.9, span = 1.2) + 
  theme_bw() + 
  xlab('\nIndependence\n(Number of Standard Deviations from Mean)\n') + 
  ylab('Patronage-Expertise Dimension\n(Number of Standard Deviations from Mean)\n') + 
  facet_grid(Weighted~Type, scales = "free")



g1.big.exp <- ggplot(selin.plot.frame.exp,  aes(x = Dimension, y = latentpatron, 
                                                  weight = Weight)) + 
  geom_point(alpha = I(0.3)) + 
  stat_smooth(color = "black", level = 0.9, span = 1.2) + 
  theme_bw() + 
  xlab('\nIndependence\n(Number of Standard Deviations from Mean)\n') + 
  ylab('Patronage Dimension\n(Number of Standard Deviations from Mean)\n') + 
  facet_grid(Weighted~Type, scales = "free_x")

g1.big.both.exp <- ggplot(selin.plot.frame.exp,  aes(x = Dimension, y = patronexpert, 
                                                       weight = Weight)) + 
  geom_point(alpha = I(0.3)) + 
  stat_smooth(color = "black", level = 0.9, span = 1.2) + 
  theme_bw() + 
  xlab('\nIndependence\n(Number of Standard Deviations from Mean)\n') + 
  ylab('Patronage-Expertise Dimension\n(Number of Standard Deviations from Mean)\n') + 
  facet_grid(Weighted~Type, scales = "free")



ggsave(g1.big, file = "Figure-D2.pdf",
       width = 7, height = 4)

ggsave(g1.big.both, file = "Figure-D5.pdf",
       width = 7, height = 4)

ggsave(g1.big.orig, file = "Figure-D1.pdf",
       width = 7, height = 4)

ggsave(g1.big.both.orig, file = "Figure-D4.pdf",
       width = 7, height = 4)

ggsave(g1.big.exp, file = "Figure-D3.pdf",
       width = 7, height = 4)

ggsave(g1.big.both.exp, file = "Figure-D6.pdf",
       width = 7, height = 4)



### pca plots
patron.dim2 <- patron.dim

row.names(patron.dim2$ind$coord) <- 
  row.names(patron.dim2$ind$cos2) <- 
  row.names(patron.dim2$ind$contrib) <- as.character(agency.obs2$FSAgency)
  

row.names(patron.dim2$var$coord) <- c("Last Job in Politics", "Campaign or Transition", "Obama Connection")
row.names(patron.dim2$var$cor) <- c("Last Job in Politics", "Campaign or Transition", "Obama Connection")
row.names(patron.dim2$var$cos2) <- c("Last Job in Politics", "Campaign or Transition", "Obama Connection")
row.names(patron.dim2$var$contrib) <- c("Last Job in Politics", "Campaign or Transition", "Obama Connection")

patron.biplot <- fviz_pca_biplot(patron.dim2, labelsize = 3, 
                                 col.ind = "grey70",
                                 col.var = "black",
                                 axes.linetype = "dotted") + 
                      xlab(paste("\nFirst Principal Component\n(", 
                                 sprintf("%2.2f", patron.dim2$eig[1,2]),
                                 "% of Variance Explained)",sep="")) + 
                      ylab(paste("Second Principal Component\n(", 
                                 sprintf("%2.2f", patron.dim2$eig[2,2]),
                                 "% of Variance Explained)\n",sep="")) +
                      xlim(-6.5, 6.5) +
                      ggtitle("")
  
ggsave(patron.biplot, file = "Figure-A1.pdf", width = 6, height = 6)

both.dim2 <- both.dim
both.dim2$ind$coord[,1] <- -both.dim2$ind$coord[,1]
both.dim2$var$coord[,1] <- -both.dim2$var$coord[,1]

row.names(both.dim2$ind$coord) <- 
  row.names(both.dim2$ind$cos2) <- 
  row.names(both.dim2$ind$contrib) <- as.character(agency.obs2$FSAgency)


row.names(both.dim2$var$coord) <- c("Postgraduate Degree", 
                                    "Bachelors Degree Only", 
                                    "Bush or Clinton Experience", 
                                    "Agency Experience", 
                                    "Subject Experience", 
                                    "Government Experience",
                                    "Last Job in Politics", 
                                    "Campaign or Transition", 
                                    "Obama Connection")
row.names(both.dim2$var$cor) <- c("Postgraduate Degree", 
                                  "Bachelors Degree Only", 
                                  "Bush or Clinton Experience", 
                                  "Agency Experience", 
                                  "Subject Experience", 
                                  "Government Experience",
                                  "Last Job in Politics", 
                                  "Campaign or Transition", 
                                  "Obama Connection")
row.names(both.dim2$var$cos2) <- c("Postgraduate Degree", 
                                   "Bachelors Degree Only", 
                                   "Bush or Clinton Experience", 
                                   "Agency Experience", 
                                   "Subject Experience", 
                                   "Government Experience",
                                   "Last Job in Politics", 
                                   "Campaign or Transition", 
                                   "Obama Connection")
row.names(both.dim2$var$contrib) <- c("Postgraduate Degree", 
                                     "Bachelors Degree Only", 
                                     "Bush or Clinton Experience", 
                                     "Agency Experience", 
                                     "Subject Experience", 
                                     "Government Experience",
                                     "Last Job in Politics", 
                                     "Campaign or Transition", 
                                     "Obama Connection")

both.biplot <- fviz_pca_biplot(both.dim2, labelsize = 3, 
                                 col.ind = "grey70",
                                 col.var = "black",
                                 axes.linetype = "dotted") + 
  xlab(paste("\nFirst Principal Component\n(", 
             sprintf("%2.2f", both.dim2$eig[1,2]),
             "% of Variance Explained)",sep="")) + 
  ylab(paste("Second Principal Component\n(", 
             sprintf("%2.2f", both.dim2$eig[2,2]),
             "% of Variance Explained)\n",sep="")) +
  xlim(-20, 15) +
  ggtitle("")

ggsave(both.biplot, file = "Figure-A2.pdf", width = 6, height = 6)




hhl.patronexpert.weight <- ols(patronexpert ~ clintonlewisideo + lnworkforcesize + 
                                                      lnprobias + 
                                                      sotupriority, 
                                        data = selin.agency.merged, x = TRUE, y = TRUE,
                                        weights = numappointees)

hhl.latentpatron.weight <- ols(latentpatron ~ clintonlewisideo + lnworkforcesize + 
                                         lnprobias + 
                                         sotupriority, 
                                       data = selin.agency.merged, x = TRUE, y = TRUE,
                                       weights = numappointees)

texreg(list(hhl.latentpatron.weight,
            hhl.patronexpert.weight),
       stars = c(0.02, 0.10, 0.20),
       digits = 3,
       custom.coef.names = c("Constant", "Agency Conservatism",
                             "Agency Workforce Size",
                             "Professionalism",
                             "Priority Agency"),
       custom.model.names = c("Patronage-Only Dimension", "Patronage-Expertise Dimension"),
       booktabs = TRUE, dcolumn  = TRUE, table = FALSE,  use.packages = FALSE,
       file = "Table-F1.tex")
